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1 Introduction 



In the de Broglie-Bohm pilot- wave formulation of quantum theory [1-10], the 
quantum state of an individual system is supplemented by a deterministic tra- 
jectory in configuration space. It is usually assumed that an ensemble of sys- 
tems with wave function ip^q^t) (in configuration space) has configurations q 
distributed according to the Born probability rule p{q,t) = \ip{q,t)\'^. More 
precisely, it is assumed that 0) — \'ip{q,0)\'^ at some initial time t = 0, 
the equations of motion then guaranteeing that this quantum distribution is 
preserved at all later times. 

The status of this assumption, that p = \'ip\'^ at some initial time, has been 
the subject of debate. Soon after Bohm's first papers on the theory appeared 
12], Pauli objected that taking a particular distribution p = \tp\'^ as an initial 
condition was unjustified in a fundamentally deterministic theory |llj . Simi- 
larly, it was pointed out by Keller 12 that p — should be derived from 
the dynamics, in order for quantum theory truly to emerge as the statistical 
mechanics of an underlying deterministic theory. 

Despite these criticisms, Bohm had in fact already suggested in his original 
papers that the distribution p — had a status similar to that of thermal 
equilibrium in ordinary statistical mechanics, and he had conjectured that this 
distribution could be derived by appropriate statistical-mechanical arguments. 
In a subsequent paper, Bohm studied the specific case of an ensemble of two- 
level molecules subjected to random external collisions, and argued that an 
initial nonequilibrium distribution p ^ would relax to the equilibrium state 
p = '13'. A general argument for relaxation was not given, however. A year 
later, in response to the criticisms of Pauli and Keller, and motivated by the 
difficulties in formulating a relaxation argument for an arbitrary system, Bohm 
and Vigier modified the theory by introducing stochastic 'fluid fluctuations' that 
drive the required process p — > \ip\'^ for a general system jl4) . 

Most subsequent presentations of the original (deterministic) theory have 
taken p — l?/;^ as an axiom. Usually, it is simply assumed that any ensem- 
ble prepared at time t = with wave function 'ip{q,0) will have a distribution 
p{q,0) — |V'(9, 0)p of actual configurations q. Alternatively, has been re- 
garded as the natural measure of probability or 'typicality' for initial config- 
urations of the whole universe (taking ijj to be the universal wave function), 
yielding the Born rule for all subsystems at all times |15| . 

However, a quite general argument for the relaxation p —> \ip\'^ may be framed 
in terms of an analogue of the classical coarse-graining 7/-theorem, based on the 
i7-function 



(equal to minus the relative entropy of p with respect to IV'P) ^IQl^j. Like 
its classical counterpart, this theorem provides a general mechanism in terms 
of which one can understand how equilibrium is approached, while not proving 
that equilibrium is actually reached (see Sect. 5 below). 




(1) 
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A relaxation timescale r - analogous to the scattering time of classical kinetic 
theory - may be defined in terms of time derivatives of the coarse-grained H- 
function. For a particle of mass m, whose wave function has a quantum energy 
spread A£', the timescale is estimated to be [T7] 



emi/2(A£;)3/2 

where e is the coarse-graining length. One expects t to be the timescale over 
which there will be a significant approach to equilibrium. This rough estimate 
agrees quite closely with numerical simulations for particles moving in one spa- 
tial dimension [T7] . 

In this paper, we study the relaxation p — > \^\'^ numerically, for the more 
realistic case of particles moving in two spatial dimensions. We shall see that 
if ?/) is a superposition of energy eigenstates, the de Broglie-Bohm velocity field 
varies extremely rapidly in general, particularly in regions where |'0| is small. 
The resulting particle trajectories have a very complicated structure. Starting 
from a simple nonequilibrium distribution po 7^ IV'oP at t = 0; the relaxation 
P ^ IV'P takes place rather efficiently, on a timescale of order And a plot of 
the coarse-grained 77-function against time shows an approximately exponential 
decay, with time constant of order 0. 

We interpret these results as further evidence that, in the pilot-wave formula- 
tion of quantum theory, the Born distribution p = should not be regarded 
as an axiom. Rather, it should indeed be seen as dynamically generated, in 
the same sense that one usually regards thermal equilibrium as arising from a 
process of relaxation based on some underlying dynamics. 

In Sect. 2, we introduce pilot-wave dynamics for a single particle in a two- 
dimensional box. In Sect. 3, we illustrate the character of the trajectories. 
In Sect. 4, we introduce a nonequilibrium ensemble of systems, and show in 
detail how the ensemble evolves towards p = IV'I^- In Sect. 5 we show how the 
relaxation may be quantified in terms of the coarse-grained _ff-function, which is 
found to decay approximately exponentially with time, on a timescale of order 
((5J. In Sect. 6, we discuss the significance of these results, addressing some 
relevant issues in the foundations of statistical mechanics. 



2 Pilot- Wave Dynamics in Two Dimensions 

Consider a single particle of unit mass, moving in a potential V in two spatial 
dimensions. The system has a configuration q — (a;,y), and wave function 
ij) = ip(^x,y,t) (assuming a pure state) satisfying the Schrodinger equation 

(units h= 1). 
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In de Broglie-Bohm theory, the particle traces out a definite trajectory 
{x{t),y(t)), which is determined by the wave function ip according to the de 
Broghe guidance law 

— = Im — = V5 4) 

dt V 

(where tp = \ijj\e'^). Mathematically, Im(V'0/V') is just the ratio of the quantum 
probability current to the quantum probability density. Physically, however, ^p 
is here interpreted as an objective field (in configuration space), guiding the 
motion of a single system. 

Indeed, so far we have said nothing about probabilities or ensembles, and the 
equations ||2J) and 10} define a deterministic dynamics for individual particles. 
Given the initial wave function ijj{x,y,0), Q determines the wave function 
iplx.yjt) at all times; and given the initial particle position (a;(0), ?/(0)), 
then determines the particle trajectory {x{t)^y{t)) at all times. 

For an ensemble of independent particles, each guided by the same wave 
function ip{x,y,t), we may define a density p(x,y,t) of actual configurations 
{x,y) at time t. The guidance equation defines a velocity field (i,y), and 
for an ensemble of particles with the same wave function the field (i, y) 
determines the evolution of any distribution p{x, y, t) via the continuity equation 

5p ^ d(px) ^ djpy) ^ ^ 
dt dx dy 

Because implies the continuity equation 

dt dx dy 

for \tp\^ , the particular initial distribution p{x,y,0) = \il!{x,y,0)\'^ evolves into 
p(x, y, t) — |^(a;, y,t)\ . The state of 'quantum equilibrium' is preserved by the 
dynamics. 

Note that Q determines the evolution of any initial distribution p{x,y,Q), 
even if p{x, y, 0) ^ li'ix, y, 0)p. Given %p{x, y, 0), (jSjl determines ^{x, y, t) at all 
times, and then determines the velocity field (i, y) at all times. Once the 
velocity field (i, if) is known everywhere, Q may be integrated to yield p(x, y, t) 
at all times, for any initial p{x, y, 0). 

Before we proceed, it is instructive to compare the above with the analo- 
gous classical Hamiltonian evolution on phase space. Classically, the trajectory 
{q{t),p{t)) of an individual system is determined by Hamilton's equations 

. _dH_ . _ dH 
dp ' ^ dq 

given the initial values (gojPo)- These equations define a velocity field {q,p) 
on phase space, and for an ensemble of systems with the same Hamiltonian 
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H the field {q,p) determines the evolution of any distribution p{q,p,t) via the 
continuity equation 

^ ^ djpq) ^ djpp) ^ ^ 
dt dq dp 

which may be rewritten as 

§ + {P,^} = (9) 

Because dq/dq = —dp/dp, a uniform initial distribution p(q,p,0) = const, (on 
the energy surface) remains uniform, p(q,p,t) — const.. The state of thermal 
equilibrium is preserved by the dynamics. While for a general (non-uniform) 
initial state p{q,p,0), the evolution p{q,p,t) is obtained (in principle) by inte- 
grating (0. 

It is known that in appropriate circumstances, the classical evolution on 
phase space defined by © leads to thermal relaxation on a coarse-grained level 
(see Sect. 5 below). Less is known about the corresponding evolution (on 
configuration space) defined by ((S)). 

Classically, the canonical example of thermal relaxation (for an isolated sys- 
tem) is that of a simple initial distribution p{q,p,0) (with no fine-grained mi- 
crostructure) concentrated on a small region of the energy surface. Under the 
Hamiltonian evolution (9), the distribution p{q,p,t) will generally develop a 
complex filamentary structure that spreads out over the energy surface, so that 
the distribution approaches uniformity on a coarse-grained level. 

Here, we shall consider an analogous example, of particles in a two-dimensional 
box with a simple initial distribution p{x,y,Q) ^ |'0(a;,y, 0)p (with no fine- 
grained microstructure). We shall see, by numerical integration of (jU, that the 
ensemble evolves towards the equilibrium distribution (on a coarse-grained 
level) . 

Consider, then, a particle confined to a square box of side tt, with infinite 
barriers at x, y — Q, tt. The energy eigenf unctions are 

2 

0mn (x, y) = - sin (mx) sin (ny) (10) 

TT 

with energy eigenvalues Emn — ^{rri^ + n^), where m,n = 1, 2, 3, are positive 

integers. 

As a specific example, we take the initial wave function ip{x,y,0) at t = 
to be a superposition of the first 16 modes, m, n = 1, 2, 3, 4, with amplitudes of 
equal modulus but randomly-chosen phases 0„ 



ip{x,y,0)^ ^ -0m«(a;,?/)exp(?6'„„) (11) 

m,n— 1 

The squared-amplitude \ip{x,y,0)\'^ is shown in Fig. 1. 

^For the record, the values of 9mn used were (to four decimals) 611 = 5.1306, 812 = 2.0056, 

6)13 = 4.1172, 014 = 3.3871, 6*21 = 6.2013, 622 = 4.6598, 623 = 1.8770, 624 = 4.3033, 

031 = 4.0145, 032 = 6.1142, 6»33 = 5.4401, 6*34 = 1.9292, 641 = 3.4015, 642 = 6.2109, 
043 = 6.0370, 6*44 = 5.9159. 
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Figure 1: The squared-amplitude \ip{x,y,0)\'^ at t = 0, for the specified super- 
position of the first 16 modes of a two-dimensional box. 



From the wave function %l:{x,y^t) at later times is just 

ip{x,y,t)= ^ -(l)^n{x,y)expi{d„in - Emnt) 



(12) 



m,n— 1 



Note that ip is periodic in time, with period An (since AirEmn is always an integer 
multiple of 2tt). 

From |@J), the velocity components of the particle at any moment are given 

by 

dx i ( , dip* I ^ dip^ 
dx 



^ dx 



and 



dt 2\tP\ 
dt 2|V'p V dy dy 



(13) 
(14) 
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3 Character of the Trajectories 



The velocities ()13|l and H14() are ill-defined at nodes (where \'4)\ = 0), and tend 
to diverge as nodes are approached. This is because, close to a node, small 
displacements in x and y can generate large changes in phase S = Imln?/;, 
corresponding to a large gradient VS".^ Note, however, that as shown by Berndl 
et al. |18[ [T^ , for reasonable Hamiltonians and wave functions the set of initial 
particle positions that reach singularities of the velocity field in finite time is of 
IV'p-ineasure zero. 

Thus, for almost all initial values (a::(0), i/(0)), the trajectory {x{t),y{t)) of 
the particle may be calculated by numerical integration of and (|14|l . 
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Figure 2: A typical trajectory, with an apparent cusp at C. The particle actually 
turns around very slowly at C. 

Our method is based on the Runge-Kutta-Fehlberg algorithm 20 , with an 
adaptive time step h adjusted so that the absolute errors in x{t) and y{t) at 

^Because ^ is a smooth, single-valued function, a small displacement (Sx, Sy) produces a 
small change (5i/) in the complex plane. However, close to a node (li/ij = 0), (SV liss near the 
origin of the complex plane and so can correspond to a large change SS of phase. 
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each step are both less than hA for some fixed A. For each trajectory, we begin 
by setting A = 10~^. We then recalculate the same trajectory with A = 10~^. 
If the global error (the distance between the final positions) is less than 0.01, 
wc retain the more accurate result. Otherwise, we recalculate with A = 10~^, 
and so on. 

For a small number of trajectories, an accurate calculation was beyond the 
bounds of our computational resources, which wc fixed at a minimum of A = 
10~^^ and a maximum of 10* time steps (per trajectory). Over a time 27r the 
number of problematic trajectories is about 1 in 60,000, while over a time iw it 
is about 1 in 1,500 (sampling uniformly over the box). 

A typical trajectory is shown in Fig. 2. In general, the trajectories are rather 
irregular. 

At some points, for example at C (Fig. 2), there is an apparent cusp, but 
closer examination shows that the tangent to the curve is not discontinuous. 
Further, at least in this case, the velocity field at C is in fact very small: the 
particle is slowly turning around, with a speed in the range ~ 10~^ — 10~^. 
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Figure 3: Close-up of a trajectory near a node or quasi-node (where is known 
to be very small and possibly zero). The particle rapidly circles around a moving 
point at which l/\tp\ is highly peaked. 
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The motion tends to be particularly rapid in regions where is small. Fig. 
3 shows a close-up of a trajectory near a nodal or quasi-nodal point, where I'^l is 
known to be very small (but not known to be strictly zero, it being impossible to 
tell in a numerical analysis). The spatial region shown is about 0.3% of the area 
of the whole box, and the displayed trajectory covers a time interval At = 0.23.^ 
The particle follows a rapid circular motion around a point moving from right 
to left in the figure and the moving point is a node or quasi-node, at which 
is highly peaked. For the case shown, l/|f/'| has an apparent peak value 
of about 60, and the speed of the particle varies widely from ~ 10 to ~ 50. 




X 



Figure 4: Illustration of the divergence of neighbouring trajectories. 

In Fig. 4, we illustrate the divergence of neighbouring trajectories. Two 
distinct but nearby initial positions at t = evolve after time t = Aw into 
widely-separated final positions. For the example shown, the initial distance 
between the points is 0.005, while the final distance is 1.403. 

The divergence of neighbouring trajectories and the behaviour of Lyapunov 
exponents in pilot-wave dynamics has been studied by a number of authors 

^Thc region shown is x € (1.35,1.6), y 6 (1.21,1.33), and the time interval is t 6 
(9.39,9.62). 
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I^ESESII^. Chaotic trajectories have been reported for a range of systems: 
the two-dimensional anisotropic harmonic oscillator the quantum kicked 

rotor 1^123], the quantum cat map [221, the hydrogen atom in an external field 
|27|. the Henon-Heiles oscillator [JS], and a two-particle stationary state |29| . 

For the two-dimensional square box, chaotic trajectories have been reported 
by Frisk [51], who emphasises the importance of nodes in generating chaotic 
motion, and whose numerical simulations suggest a proportionality between the 
Lyapunov exponent and the number of nodes. 



4 Relaxation to Quantum Equilibrium 

We now turn to the time evolution of an initial nonequilibrium ensemble. We 
assume that every particle in the ensemble is guided by the same wave function 
tl>{x,y,t), given by (|12|l . At t = 0, the (nonequilibrium) distribution is chosen 
to be 



2X2 



p{x,y,0) — { — ] sin xsin y (15) 



This is just \(j)ii{x,y)\'^ - that is, the ground-state equilibrium distribution. It 
is displayed in Fig. 5. 

Instead of integrating the continuity equation (jSJ directly, it is more conve- 
nient to calculate the trajectories, from which it is straightforward to deduce 
the evolution of any initial distribution p(x, y, 0). 

From (0 and it follows that the ratio 



is conserved along trajectories: df /dt — ot 

dj_ ^ dxdf_ _^ dydf_ ^ ^ 
dt dt dx dt dy 

If the initial point (a;o, yo) evolves into (x, y) at time then for any ensemble 

/(x,y,i) = /(xo, 2/0,0) (17) 

Thus, given the function 

f[XQ,yQ,0) = T—f ?VU2 

|i/;(a;o,yo,0)|^ 

at t = 0, the distribution p(x, y, t) at later times may be written as 

p{x,y,t) = \ip{x,y,t)ff{xo,yo,0) (18) 

Given the mapping from (xq, yo) to {x,y), the fate of any initial distribution 
p{x,y,0) may be deduced from (fl^ . 
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Figure 5: The initial nonequilibrium distribution, chosen to be p{x,y,0) = 
|(/>ii(x, ?/)|^ (equal to the ground-state equilibrium distribution). 



One might begin with a uniform lattice of initial points {xo,yo), and cal- 
culate the future trajectories {x{t),y{t)), using H18|) to deduce the distribution 
p{x{t),y{t),t) at the later points {x{t),y{t)). This method was used in the first 
relaxation simulation, for the simple case of a one-dimensional box However, 
the lattice of points at later times is unfortunately distorted, and there appear 
large regions containing hardly any lattice points at all. 

A more accurate procedure is to evolve trajectories backwards, from (x, y) 
at any desired t to the starting point (a;o,?/o) at i = j30|. Thus, to calculate 
p{x,y,t) at time t, we set up a uniform lattice of points {x,y) at time t and 
backtrack the trajectories to the initial points {xo,yo) at i = 0. Knowing the 
function f{xo, yo, 0), we may then use H18|l to deduce y, t). 

This backtracking procedure has the disadvantage that one has to calculate 
the trajectories over the whole interval from t to t = 0, for every desired value of 
t (for which one wishes to know p(x, y, t)). But overall the procedure is better, 
because it generates results for p{x, y, t) on a uniform lattice. This is particularly 
important when one wishes to calculate the coarse-grained 7J-function. 
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As mentioned in Sect. 3, for some trajectories an accurate calculation is 
beyond our computational resources (A = 10~^^ and 10^ time steps per trajec- 
tory). Nevertheless, for practical reasons of data handling, it is convenient to 
assign values to all of them, rather than ignore the problematic cases. ^ Since 
the number of problematic trajectories is found to be small (1 in 60,000 over 
a time 27r and 1 in 1,500 over a time Att), they have a negligible effect on the 
results reported below. 

Note that because the values of / = p/\tp\'^ are carried along trajectories, 
the exact (fine-grained) value of p{x,y,t) will always differ from \ilj{x,y,t)\'^, 
by just the same multiplicative factor f{xo,yo,0) by which p(.To,j/o,0) differed 
from \ip{xo, yo) 0)p. Indeed, the property df /dt = is analogous to the Liouville 
property in classical statistical mechanics, according to which the phase-space 
density is constant along Hamiltonian trajectories (sec Sect. 5 below). In both 
cases, relaxation occurs only in a coarse-grained sense. 

To see how this works, in Fig. 6 we show a close-up of p{x, y, t) in a small 
fixed region of the box, at time t = n. The exact, fine-grained distribution is 
highly irregular. The area shown is very small: it has side 7r/128, covering just 
0.006% of the total area of the box. (For this calculation, we used a uniform 
lattice of 40, 000 points in this small area, with the lattice specified at t = tt.) 

Clearly, the exact fine-grained distribution is extremely complex, and it is 
useful to consider its appearance under coarse-graining. Given a (square) coarse- 
graining cell of side e, the distribution p may be averaged over the cell to produce 
a coarse-grained value 



which may be assigned to the centre of the cell. 

If the cells do not overlap, we denote the coarse-grained distribution by p. 
Results for p will be shown in Sect. 5 below (using a coarse-graining length 
s = 7r/32). 

A smoother coarse-grained distribution, denoted p, is obtained if one takes 
overlapping cells (a procedure used in quantum chemistry, for example). To 
generate p, we take cells of side e = n/ld that overlap with their immediate 
neighbours in 88% of their area. Specifically, given one cell, shifting it along x 
or y by a distance equal to 12% of its length (that is, 12% of 7r/16) generates 
the neighbouring cell, and so on. In this way, in the box of area tt^ we obtain 
a patchwork of 126 x 126 overlapping cells, each with its own value of p. For 
the trajectory calculations, we use a lattice of 400 x 400 particles, with 25 x 25 
particles in each cell. 

We have calculated the evolution of p up to t = Aw (when recurs to 
its initial value). This is shown in Fig. 7. Displayed are p and at times 

*Some cases required more than 10* time steps even for A = 10~®; each of these were 
aborted and assigned a value from a neighbouring trajectory. In other cases, decreasing A 

exceeded tiic limit of 10* time steps, and we retained the value obtained using the smallest 
value of A that did not exceed this limit (despite the global error exceeding 0.01). Finally, 
some cases reached the minimum setting A = 10^ ^'^ without achieving the desired limit on 
the global error, and for these we assigned the value obtained using A = 10"^^. 




cell 



dxdy p 



(19) 



12 



0.7 n 



o.e- 




y 



Figure 6: Close-up of the exact p{x^ y, t) in a small spatial region, covering just 
0.006% of the area of the box, at time t = tt. The fine-grained distribution is 
very irregular. 



i = 0, 27r and An. Evidently, p and become rather close, and a remarkable 
relaxation towards quantum equilibrium has occurred after just one time cycle 
of 47r. 

Fig. 8 shows the same data, as contour plots. 

Note that while ip and its associated velocity field are periodic, in general 
the trajectories are not, so that the distribution p does not recur at t = in. 
This is in contrast with the one-dimensional case, where the particles cannot 
move past each other, constraining each trajectory to recur simultaneously with 
IV'P, so that any initial po ^ IV'oP recurs as well piTT]. 

Presumably, evolution over more time cycles of Aw will make p converge ever 
closer to |V'p. However, we have not extended our calculations beyond t = An. 
Possibly, a small residual nonequilibrium will remain after an arbitrary number 
of cycles, but this seems unlikely. 

Given the mapping from {xo,yo) to (a;(47r), ?/(47r)), not only can one im- 
mediately construct p{x,y^AT:) from any initial p{x,y,Q)^ one can also extend 
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Figure 7: Smoothed p, compared with at times t = 0, 27r and 47r. While 
I^P recurs to its initial value, the smoothed p shows a remarkable evolution 
towards equilibrium. 



the evolution to an arbitrary number of cycles An that is, to times t = Ann 
for positive integral n - by simple iteration of the mapping (the velocity field 
being periodic with period 47r). It might then appear that, given the calcula- 
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Figure 8: Smoothed p, compared with \ip\'^, at times t = 0, 2n and Air. The 
same data as in Fig. 8, displayed as contour plots. 



tions performed so far, we could easily extend them to t = 47rn. Unfortunately, 

the calculated map from {xqjJJq) to (a;(47r), i/(47r)) cannot be reapplied to map 
(x(47r), j/(47r)) to (x(87r), 2/(877)) and so on, because the backtracking method we 
have used distorts the lattice at t = 0. 
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5 Decay of the iif-Function; Relaxation Timescale 

The difference between p and may be quantified by the value of the H- 
function As we mentioned in the Introduction, this is just minus the relative 
entropy of p with respect to jV'l^ (a standard measure of the difference between 
two distributions). The 77-function is bounded below by zero, and is equal to 
zero if and only if p — | | ^ everywhere TB", T, 17' . 

Note, however, that because of (|SJ| and © the exact H is constant in time, 
dH/dt = 0, reflecting the fact (already mentioned) that differences between fine- 
grained values of p and never disappear (the ratio / = p/\'4>\'^ being carried 
along trajectories). A decreasing H is obtained only under coarse-graining. 

The situation is similar classically, where departures from thermal equilib- 
rium may be quantified in terms of the classical _ff-function (on phase space) 



which is just minus the relative entropy of p with respect to the uniform distri- 
bution. Under Hamiltonian evolution, dp/dt — along trajectories (Liouville's 
theorem), and the exact ffciass is constant, dHdass/dt = 0. However, in appro- 
priate circumstances, the classical phase-space evolution defined by Q leads to 
thermal relaxation on a coarse-grained level. This may be quantified in terms 
of the classical coarse-grained 77-function 



where p is p averaged over (non-overlapping) coarse-graining cells. This obeys 
the coarse-graining 77-theorem |^ 15^ 



(assuming no initial fine-grained microstructure in p at t = 0), and i?ciass is min- 
imised by p = const.. This theorem formalises the intuitive idea of Gibbs - that 
an initial non-uniform distribution will tend to develop fine-grained structure 
and become more uniform on a coarse-grained level. 

In the pilot-wave case too, the evolution (on configuration space) may lead 
to relaxation on a coarse-grained level, as shown by the above numerical simu- 
lations. And the relaxation may be similarly quantified in terms of the coarse- 
grained if-function 





-^class (i) < 5elass(0) 




(20) 



which obeys the _ff-theorem [161 0] El 



H{t) < H{0) 



(assuming no initial fine-grained microstructure in p and IV'I^); and where again 
H > for all p, and H = Q ii and only if p = \tp\^ everywhere. This 
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version of the H -theorem formahses the idea that p and \tjj\^ behave Hke two 
'fluids' which are 'stirred' by the same velocity field (since they obey the same 
continuity equation), so that p and tend to become indistinguishable on a 
coarse-grained level. ^ 

A natural timescale r for relaxation may be defined in terms of the time 
derivatives of H{t) at t = 0. 

Using the continuity equations |SJ| and ©, it follows that (^dH/dt)^ — and 
(for g= (a;,2/)) 01 



|2 



dxdy^ {Xo ■ V/o)2 - {Xo ■ V/o) < (21) 



dt^ Jo J J fo 

where = 'tp{x,y,0), fo = f{x,y,0), Xo = {io,yo) and V = {d/dx,d/dy). 
(The quantity in brackets is non-negative, and is equal to the spatial variance 
of Xo ■ V/o over a coarse-graining cell.) 
Thus we may define r by 01 

where H21f) gives {cPH /dt^^^ in terms of the initial conditions ipo and po- 

Expanding Xo ■ V/o in a Taylor series within each coarse-graining cell of area 
e^, it is found that 

£2 



[d^H/de)^^~-I + 0{e^) (23) 



where [HI [101 



Thus 



%^|V(Xo- V/o)|' (24) 
Jo 



r-\\l'^ + 0{e) (25) 

If £ is small with respect to the lengthscale over which Xq ■ V/o varies, then 
T oc l/£. Taking Ho 1 (a mild nonequilibrium) , a rough estimate of / yields 
the quoted result in terms of the quantum energy spread Ai5 of ■00 [TTIITnj . 

Note that, given the result r oc l/£ (for small £), the formula ^ may be 
obtained on dimensional grounds. And in any case, (O is only a rough, order- 
of-magnitude estimate. For our initial wave function (|ll(l . Ai? ~ 4 and (in our 
units) 

.^1 (26) 

Thus, at least initially, H — H{t) should decay on a timescale of order l/8£. 

^Note, however, that the velocity field jlj is related to ij) by the Schrodinger equation, 
as well as by the continuity equation, whereas no such additional relation exists between the 
velocity field and p. This explains why p develops a filamentary structure while does not 

in 
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To investigate the decay of H{t) numerically, we must work in terms of the 

coarse-grained quantities p and {ipl^ , with non-overlapping coarse-graining cells 
(as used in the above if-theorem). 





Figure 9: Coarse-grained p and coarse-grained \tp\ (with non-overlapping cells) 
at times t — and t — 27r. 



Fig. 9 displays p and |-0| , with non-overlapping cells of side e = 7r/32, at 
times t = and t = 2-k. 

In Fig. 10, we plot In H against time from t = Qiot — 211 (where the value of 
H was calculated every 7r/4 time units). The evolution of In^ is approximately 
linear, corresponding to an approximately exponential decay H(t) « ^ffpG"*^'"! 
with a (best-fitting) time constant tc ~ 4. This compares favourably with the 
order-of-magnitude estimate ]2<j\i . which for the present coarse-graining length 
e = n/32 gives r ~ 1.3. 
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1 2 3 4 5 6 

t 

Figure 10: Plot of \nH against time, from t — Q to t — 2tt, showing an approx- 
imately exponential decay of H{t) with a time constant tc ~ 4. 



Fig. 10 shows error bars for the displayed values of InH. These were ob- 
tained as follows. The main source of error is the approximation of p by an 
average over a finite number of points within a coarse-graining cell. To estimate 
this error, each calculation of p is done again with a different lattice of points, 
corresponding to a different sampling within the cells (first with a lattice of 
25 X 25 particles, then with 27 x 27). The two results yield different values for 
H; however, the difference is usually only about 2%, and never more than 4%. 
Another (considerably smaller) source of error in H is of course the inaccuracy 
in the calculated particle positions. This error is estimated by comparing values 
of H obtained with local errors A (in the trajectories) differing by a factor of 
ten; the difference in values of H is found to be quite negligible. 

Note that the natural quantum timescale for this system is At ~ h/AE ~ 
0.25, which is arguably not much further from tc than is r. However, generally 
speaking, the quantum timescale At ~ h/AE for a system is unrelated to the 
estimated relaxation timescale t ~ {1/ e)h'^ /m}^'^ (AE)^/^ in (O, since At has 
no dependence on e and scales differently with AE. Writing AE ~ {Ap)^/2m, 
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the special choice of e S/Ap for the coarse-graining length leads to r ~ At, 
and our estimated r is not far from Ai only because our e is not far from ft/Ap. 

It would be interesting to check numerically the dependence of the relaxation 
timescale on the coarse-graining length e and on the energy spread Ai?. This has 
been done for the rather artificial but much simpler case of the one-dimensional 
box (in which trajectories cannot move past each other), with results that accord 
quite well with the predicted scahngs r oc 1/e and r oc l/(Aii;)^/^ QTlEll. To 
perform the corresponding checks for the two-dimensional case is a task for the 
future. 

We have found an approximately exponential decay of li with time. Pre- 
sumably this behaviour could be derived, by supplementing the underlying dy- 
namics with some sort of phenomenological Markovian assumption, analogous 
to the classical hypothesis of molecular chaos at every instant (from which the 
Boltzmann equation may be derived). 

6 Discussion and Conclusion 

We have seen that the simple choice of an initial nonequilibrium state (|15|l 
relaxes towards equilibrium rather efficiently, with an accompanying (approxi- 
mately exponential) decrease in the coarse-grained ff-function on the expected 
timescale. 

Of course, this is not to imply that any initial distribution po 7^ IV'oP would 
relax to equilibrium: that is impossible in any time-reversal invariant theory. 
Just as in classical mechanics, for every solution of the dynamical laws evolv- 
ing towards equilibrium, there is a time-reversed solution evolving away from 
equilibrium. For example, as new initial conditions -00, at t = we could 
take 

V''(a;,2/,0) = ?A*(x,y,7r), p'(x, y, 0) = p(a;, tt) (27) 

where p(x^ y, tt) is the former (exact) distribution at t = n, whose complex fine- 
grained microstructure is shown in close-up in Fig. 6. The dynamical equations 
((31 and © imply that these conditions will evolve into 

i''{x,y,t) = 'ilj*{x,y,7r -t), p' {x,y,t) ^ p(x,y,7r - t) (28) 

so that after a time t = n we recover our former initial distribution p{x, y, 0) 
given by (|15|l - which is further away from equilibrium than the present initial 
state H27|l . Thus the initial conditions 1)27(1 evolve away from equilibrium, and 
the coarse-grained iJ-function increases with time. 

However, as in ordinary statistical mechanics, relaxation towards equilibrium 
will be obtained if one assumes that the initial state satisfies certain conditions 
- such as the absence of fine-grained microstructure - which are violated by 
time- reversed 'initial' states such as (|27ll . Conceptually, the situation here is 
the same as in ordinary statistical mechanics. 

Conceptual issues in statistical mechanics arguably reduce to issues about 
initial conditions, and these are inevitably bound up with questions of cosmology 
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|32l I34L 051 1361 117j . According to our current understanding, we see thermal 
nonequilibrium in our universe today only because gravity has the remarkable 
property of amplifying small inhomogeneities in temperature and energy density, 
leading to the formation of large-scale structure out of a primordial homogeneous 
state of thermal equilibrium [37] . In the absence of gravity, a classical isolated 
system with no initial fine-grained structure would be expected to evolve towards 
thermal equilibrium, just as in pilot-wave dynamics. 

Initial conditions can be questioned. A lot of current work in cosmology 
is motivated by the desire to explain, for example, why the early universe was 
so smooth and homogeneous, and the prevailing view is that the homogeneity 
arose from interactions taking place at even earlier times (within cosmological 
horizons). Similarly, in the de Broglie-Bohm formulation of quantum theory, one 
may ask why all physical systems probed so far follow the Born rule p = 
Because these systems all have a long and violent astrophysical history, it is 
reasonable to explain the distribution p = ['(/'[^ in terms of a relaxation process 
from an earlier nonequilibrium state p ^ \'4'\'^ ■ This leads naturally to the 
suggestion that quantum nonequilibrium may have existed in the early universe, 
in which case the quantum noise we see today may be regarded as a remnant of 
the big bang 38 1II5 ^I17[I1(J| . A similar view may be taken in any deterministic 
hidden- variables theory PUIHTIII^ . 

In the original pilot-wave version of the i/-theorem ^S], it was argued that 
relaxation P j^fp would occur (on a coarse-grained level) for an ensemble 
of many-body systems with distribution P(xi,X2, ....X7v,t) and wave function 
\['(xi,X2, ....xjv,i) (with N large). It was then shown that single particles ex- 
tracted from the (eventual) equilibrium ensemble and prepared with wave func- 
tion ip would have a distribution p — However, it is now clear from the 
above simulations that a large number of degrees of freedom are not needed 
for efhcient relaxation to occur. Even for a single particle, relaxation will occur 
rapidly if its wave function is a superposition of even a modest number of energy 
eigenfunctions, as our numerical example shows. 

The simulations performed in this paper add substance to the view, already 
mentioned, that relaxation occurs because p and \iJj\^ evolve like two 'fluids' 
which are 'stirred' by the same velocity field. The most efhcient mixing is 
found to occur in the neighbourhood of nodes or quasi-nodes, where |^| is small. 
These points move around inside the box, rather like 'electric mixers' (or stirring 
devices) moving through a fluid, generating an efficient relaxation everywhere. 

We mentioned in Sect. 2 that the importance of nodes in generating chaotic 
motion was noted by Frisk 24j. Frisk also suggested that nodes would be 
important if the motion is to have the appropriate mixing properties required for 
complete relaxation to quantum equilibrium. The above simulations certainly 
bear out the expectation that nodes provide a particularly efficient mixing. But 
whether or not relaxation would take place even in their absence (perhaps on 
much longer timescales) is not known. 

Many details remain to be understood. We are far from a full understanding 
of pilot-wave theory as a dynamical system. General properties of the flow, such 
as ergodicity and mixing (in a rigorous sense), remain to be investigated. Still, 
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the results already in hand strongly suggest that, in de Broglie-Bohm theory, 
the distribution p = ji/jp is an equilibrium state quite analogous to thermal 
equilibrium in classical physics. 

It should be noted that the status of the Born rule has been a contentious 
issue in quantum theory generally, perhaps most notably in the many-worlds 
formulation of Everett 03 ESI- Some recent authors 0^103 base their jus- 
tification of the Born rule on Gleason's theorem 0Hj, which states that the 
Born rule is the unique probability assignment satisfying 'noncontextuality' 
- the condition that the probability for an observable should not depend on 
which other (commuting) observables are simultaneously measured. However, 
as pointed out by Bell 0H| j Gleason's noncontextuality condition is very strong, 
as it amounts to assuming that mutually-incompatible experimental arrange- 
ments yield the same statistics for the observable in question. Saunders [50] 
gives an 'operational' derivation of the Born rule (in which it is assumed that 
probabilities are determined by the quantum state alone); while Zurek [^J ap- 
peals to 'environment-assisted invariance'. Other recent derivations of the Born 
rule arise from novel axioms for quantum theory |52[ I53| . 

Like Euclid's axiom of parallels in geometry, the Born rule seems to stand 
apart from the other axioms of quantum theory, and there have been a number 
of attempts to derive it either from the other axioms or from something else. 
We have argued in this paper that, in the de Broglie-Bohm formulation of quan- 
tum theory, the Born rule has a status similar to that of thermal equilibrium in 
ordinary statistical mechanics, and should not be regarded as an axiom at all. 
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